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ABSTRACT 


This  thesis  reports  on  three  procedures  and  the 
associated  numerical  results  for  obtaining  semiconductor 
majority  carrier  concentrations  when  subjected  to  a 
temperature  sweep.  The  capability  of  predicting  the 
exhaustion  regime  boundaries  of  a  semiconductor  is  critical 
in  understanding  and  exploiting  the  full  potential  of  the 
modern  integrated  circuit.  An  efficient  and  reliable 
method  is  needed  to  accomplish  this  task.  Silvaco 
International's  semiconductor  simulation  software  was  used 
to  predict  temperature  dependent  majority  carrier 
concentration  for  a  semiconductor  cell.  Comparisons  with 
analytical  and  numerical  MATLAB-based  schemes  were  made. 
This  was  done  for  both  Silicon  and  GaAs  materials. 
Conditions  of  the  simulations  demonstrated  effect  known  as 
Bandgap  Narrowing. 
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EXCEUTIVE  SUMMARY 


This  thesis  reports  on  three  procedures  and  the  asso¬ 
ciated  numerical  results  for  obtaining  semiconductor 
majority  carrier  concentrations  when  subjected  to  a  tempe¬ 
rature  sweep.  The  capability  of  predicting  the  exhaustion 
regime  boundaries  of  a  semiconductor  is  critical  in  under¬ 
standing  and  exploiting  the  full  potential  of  the  modern 
integrated  circuit.  This  includes,  but  not  limited  to,  ICs 
designed  for  military  hardware  as  well  as  specialized 
extrinsic  semiconductors  designed  for  space-based  appli¬ 
cations.  An  efficient  and  reliable  method  is  needed  to 
accomplish  the  task  of  predicting  exhaustion  regime 
boundaries.  Silvaco  International's  semiconductor  simula¬ 
tion  software  was  used  to  predict  temperature  dependent 
majority  carrier  concentration  for  a  semiconductor  cell. 
Comparisons  with  analytical  and  numerical  MATLAB-based 
schemes  were  made.  Conditions  of  the  simulations  demon¬ 
strated  the  effect  known  as  Bandgap  Narrowing. 

The  goals  of  this  project  are  twofold.  The  first  goal 
is  to  design  a  program  procedure  to  temperature  sweep  a 
uniform  single-doped  semiconductor  with  a  two-dimensional 
Silvaco  simulation  software [1].  There  is  a  significant 
challenge  to  overcome  in  meeting  this  goal.  Specifically, 
this  is  due  to  a  limitation  in  Silvaco  software  in  that 
temperature  sweeps  cannot  be  performed  within  a  single 
"run".  The  second  goal  was  to  test  the  results  obtained  in 
Silvaco  against  one-dimensional  analytical  and  numerical 
models  implemented  in  MATLAB. 

Figure  ES  contains  six  plots  of  logio  ( Po )  versus 
1000x(l/T)  for  various  conditions  (with  and  without  band- 


xv 


gap  narrowing  (BGN)  and  from  various  sources  (numerical, 
analytic,  and  Silvaco)  for  the  case  of  doped  Silicon.  The 
model  for  band-gap  narrowing  [6]  used  here  is: 


E  (T)  =  E 

gV  J  go 


aT 2 

(T  +  J3)  ‘ 


(1) 


With  and  Without  Bandgap  Narrowing  Plot 


Figure  ES.l  Composite  Log  Plot  of  Three  Regimes  of 

Silicon . 

The  number  codes  identify  the  various  curves  as  per 
description  in  the  legend  of  the  figure.  In  Table  1 
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recommended  values  for  silicon  including  BGN  parameters  a 
and  (3  are  provided.  Note,  as  expected,  curves  from  the 
numerical  and  analytic  methods  are  indistinguishable,  for 
example,  MATLAB-based  analytic  curve  #4 (w/BGN)  overlays 
MATLAB-based  numerical  curve  #6 (w/BGN).  Furthermore,  as 
expected  the  impact  of  BGN  is  apparent  in  the  curves  at  the 
higher  temperatures.  The  difference  between  Silvaco  pre¬ 
dictions  and  the  1-dimensional  modeling  efforts  are  slight 
but  noted  for  future  investigation.  Similar  analysis  was 
completed  for  GaAs,  with  slight  variations  due  to  the 
different  material  structure  and  properties.  The  results 
for  GaAs  curves  were  similar. 


Parameter 

Value 

Ea 

0.045  eV 

Na 

10r/  #/ cm3 

Ego 

1.11  eV  (Silicon) 

9a 

4 

BGN  a 

0 . 00047 

BGN  (3 

636.0 

Table  ES . 1 :  Silicon  Parameters  used  in  Simulations 


TonyPlot  was  used  to  create  a  slideshow  movie,  to  aid 
in  visualization  of  the  process.  The  structure  files  from 
the  DeckBuild  runs  were  saved  and  TonyPlot  was  used  to 
display  the  hole/electron  concentration  of  the  lattice  at 
temperature  steps.  This  data  was  then  exported  into 
Microsoft  Office  PowerPoint. 
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I. 


INTRODUCTION 


A.  CHAPTER  INTRODUCTION 

This  first  chapter  will  provide  an  introduction  to  the 
work  accomplished  in  this  thesis,  a  literature  review,  and 
a  structural  overview  of  the  thesis. 

B .  PURPOSE 

The  purpose  of  this  research  is  to  compare  the 
performance  of  commercial  semiconductor  two-dimensional 
simulation  tools  with  one-dimensional  analytical  and 
numerical  models  for  freeze-out  and  exhaustion  under  a 
variety  of  nontrivial  conditions. 

C.  INTRODUCTION  TO  GENERAL  TOPIC 

Integrated  circuits  (ICs)  are  specified  to  operate 
between  designated  temperature  limits.  The  circuit  designer 
selects  the  doping  level  or  levels  and  typically  assumes 
that  the  dopants  are  approximately  100%  ionized,  i.e., 
exhaustion  and  temperature  are  not  too  high.  There  can  be  a 
significant  impact  on  a  plethora  of  device  parameters  such 
as  depletions  widths  and/or  FET  threshold  voltages  if  the 
assumption  is  violated.  In  the  domain  where  the  tempera¬ 
ture  is  too  low,  the  percentage  ionization  of  dopant  or 
dopants  will  be  significantly  less  than  100%.  The  value 
for  the  majority  carrier  concentration  is  depressed 
significantly  below  the  design  value.  On  the  other  hand,  if 
the  temperature  is  too  high,  the  thermal  generation  effect 
causes  the  majority  carrier  concentration  to  become 
excessively  higher  than  the  design  value  in  what  is  called 
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the  intrinsic  temperature  regime.  The  exhaustion  regime 
lies  between  this  two  extremes,  intrinsic  and  freeze-out. 

It  is  well  known  that  for  a  multiple  impurity  dopant 
process  this  "simple"  three-regime  description  can  be 
inadequate.  What  is  important  to  the  designer  is  the  plot 
of  the  majority  carrier  versus  the  temperature  or,  what  is 
more  commonly  done,  a  plot  of  majority  carrier  concen¬ 
tration  versus  the  reciprocal  of  temperature. 

D .  GOALS 

The  goals  of  this  thesis  are  twofold.  The  first  goal 
was  to  design  a  program  procedure  to  temperature  sweep  a 
uniform  single-doped  semiconductor  using  a  two-dimensional 
Silvaco  simulation  software [1].  There  is  a  significant 
challenge  to  overcome  in  meeting  this  goal.  Specifically, 
there  is  a  limitation  in  Silvaco  software  in  that  temp¬ 
erature  sweeps  cannot  be  performed  within  a  single  "run". 
The  second  goal  was  to  test  the  results  obtained  in  Silvaco 
against  one-dimensional  analytical  and  numerical  models 
implemented  in  MATLAB. 

E.  THESIS  OVERVIEW  AND  LITERATURE  REVIEW 

This  research  begins  in  Chapter  II  with  a  review  of 

the  analysis  of  the  different  portions  of  the  project.  The 

analysis  in  Chapter  II  serves  as  the  background  for  the 

MATLAB  numerical  and  analytic  methods.  Chapter  III  is  an 

introduction  and  description  for  the  Silvaco-based 

computational  tools  used  to  evaluate  and  compare  with  the 

MATLAB-based  methods  briefly  covered  in  Chapter  II.  The 

procedure  for  modeling  of  semiconductors  in  Silvaco  and 

other  technical  aspects  will  be  discussed.  Chapter  IV  is 
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the  presentation  and  discussion  of  the  results  obtained 
from  simulations  and  testing.  This  data  will  be  used  to 
crosscheck  the  three  approaches  to  obtaining  temperature 
dependence  of  majority  carrier  concentration  in  singly 
doped  semiconductors.  Chapter  V  is  a  submittal  of 
conclusions  on  the  study  completed  and  recommendations  for 
further  work. 

Several  topics,  which  are  discussed  in  the  appendices, 
have  been  included  for  completeness.  The  corresponding 
information  is  presented  this  way,  in  order  to  avoid 
cluttering  and  overcrowding  of  the  main  text.  Appendix  A 
provides  a  brief  introduction  to  general  semiconductor 
principles.  Appendix  B  provides  an  operational  guide  to 
reproducing  the  results  in  Silvaco,  along  with  an 
abbreviated  version  of  the  Silvaco  DeckBuild  source  code 
created  for  this  work.  Appendix  C  includes  a  pre-print  of 
the  paper  presented  at  the  45th  IEEE  International  Midwest 
Symposium  on  Circuits  and  Systems  in  Tulsa,  OK,  on  Aug  4-7, 
2002  . 

Currently,  commercial  software  packages,  such  as 
available  through  Silvaco  International  [1],  are  well 
designed  to  solve  the  electron/hole  transport  problem. 

This  type  of  calculation  is  usually  required  to  predict  the 
device  I-V  (current-voltage)  characterist ics .  Surprisingly, 
using  the  same  package  to  obtain  a  temperature  dependent 
plot  for  majority  carrier  concentration  for  a  uniform 
semiconductor  requires  a  somewhat  complicated  procedure 
[2] .  This  thesis  presents  an  efficient  and  novel  way  of 
obtaining  this  curve  from  the  Silvaco  International 
software  and  compares  the  results  with  a  proposed  one- 


3 


dimensional  single-equation  analytic  model  and  a  one¬ 
dimensional  numerical  model  that  predicts  the  temperature 
dependence  for  majority  concentration  in  all  regimes  [3, 

4] .  To-date,  most  "analytic"  methods  for  determining 
dominant  features  in  such  plots  make  use  of  multiple  ad-hoc 
arguments,  which  taken  one  at  a  time,  applies  in  only  two 
of  the  three  regimes  mentioned  above  [5,  6,  7] .  Numerical 
methods  based  on  one-dimensional  analysis  provide 
significant  flexibility  in  terms  of  making  predictions  when 
there  are  multiple  dopants  and  when  taking  into  account 
second  order  effects  [4] .  It  can  be  argued  that  numerical 
methods  are  generally  going  to  be  applicable  over  a  wider 
range  of  problems  than  what  can  be  solved  with  exact 
analysis.  Two  and  three-dimensional  numerical  finite 
element  methods  provided  in  CAD  based  commercial  packages 
have  been  widely  employed  by  the  engineers  in  the 
semiconductor  processing  industry.  The  one-dimensional 
results,  both  numerical  and  analytical,  are  compared  with 
predictions  from  such  a  package  of  tools  from  Silvaco 
International  [1].  The  numerical  one-dimensional  method  is 
based  on  solving  for  zeros  in  the  charge  neutrality 
condition  [5,6,7]  iteratively  for  selected  temperatures. 

The  numerical  algorithm  involves  a  root  finding  scheme 
known  as  the  Method  of  Interval  Bisection  [8] .  The  one¬ 
dimensional  MATLAB-based  numerical  method  can  be  applied  on 
conditions  not  specifically  studied  in  this  thesis  such  as 
multiple  dopants.  Details  on  the  one-dimensional  numerical 
algorithm  are  not  covered  in  this  thesis  and  interested 
readers  are  encouraged  to  see  References  [3]  and  [4] .  The 
analytic  expression  for  the  majority  carrier  temperature 
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dependence  is  based  on  the  well-known  solutions  to  a  cubic 
equation  [ 9 ] . 

Band  gap  narrowing  [6]  is  a  well-known  effect  in  which 
the  effective  band  gap  of  the  semiconductor  decreases  with 
increasing  temperature.  Cases  with  and  without  the  bandgap 
narrowing  effect  were  simulated  to  exhibit  the  flexibility 
and  consistency  of  the  various  approaches  taken.  Many  of 
the  salient  features  in  this  thesis  have  recently  been 
summarized  in  a  short  paper  [10] .  Background  information, 
which  appears  in  Appendix  A,  was  referenced  from  these 
sources :  [11-15] . 

F.  BENEFITS  OF  STUDY 

With  the  proliferation  of  semiconductor-based  app¬ 
lications  in  modern  society,  it  has  become  critical  to 
understand  the  effects  of  temperature  on  semiconductors 
properties.  Understanding  the  nature  of  temperature 
sensitivity  before  production  and  implementation  of  a 
device  can  reduce  cost  and  risk  factors  for  a  final 
product.  For  example,  the  issue  is  particularly  evident  in 
the  case  of  space-based  applications,  due  to  the  harshness 
of  the  environment. 

G.  CHAPTER  SUMMARY 

This  Chapter  described  the  general  area  of 
investigation,  the  goals  of  the  thesis,  and  provided  a 
literature  review.  In  the  next  chapter,  the  mathematical 
formulization,  which  serves  as  a  basis  for  one-dimensional 
analytical  and  one-dimensional  numerical  scheme  for  solving 
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for  temperature-dependent  majo 
will  be  addressed. 


rity  carrier  concentrat ion. 
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II.  TECHNICAL  ANALYSIS 


A.  CHAPTER  INTRODUCTION 

This  chapter  provides  the  background  analysis  for 
predictions  of  temperature  dependent  majority  carrier 
concentrat ion  for  both  numerical  and  analytic  MATLAB-based 
approaches  (see  Figure  2-1). 


Figure  2-1:  Layout  of  Project. 


B.  BACKGROUND  ANALYSIS 

The  following  assumptions  were  made:  1)  The  semi¬ 
conductor  is  spatially  uniform,  i.e.,  no  spatial  variation 
of  the  doping  and  energy  bandgap.  2)  The  semiconductor  is 
in  equilibrium.  3)  "Degeneracy  condition"  [6],  which  is  a 
byproduct  of  excessive  doping,  has  not  occurred.  These 
assumptions  were  used  throughout  the  simulation  for  both 
MATLAB  and  Silvaco  software. 

The  Fermi-Dirac  occupation  probability  for  electrons 
and  holes  respectively  are: 
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and 


FD(E)  = 


1 

(e-Ef)  /  kT 


(2.1) 


1  +  e 


1 


FD(E)  =  1  -  FD(E)  = 


1  +  e<.EF-E)/kT 


(2.2) 


where  E  is  the  energy  level  defined  on  electron  energy  band 
diagram  (see  Figure  2-2) .  The  Fermi  energy  level,  EFj  is 
where  electron  and  hole  occupation  probability  are  both 
50%,  k  is  Boltzman' s  constant,  and  T  is  temperature  in 
Kelvin.  Figure  2-2  is  a  representation  of  the  energy 
bandgap  diagram  of  a  semiconductor.  On  the  right  hand  side 
are  the  levels  for  the  conduction  and  valance  levels,  Ec  and 
Ev.  The  total  bandgap  size  is  represented  as  Eg,  and  AEd 
and  AEa  are  used  to  represent  the  difference  between  the 
donor  (Nd)  and  acceptor  (Na)  levels  and  the  band  limits. 


Na 


Figure  2-2.  Diagram  of  Bandgap  Layout. 
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The  three-dimensional  electron  and  hole  density-of 


states  [6]  are  given  by  respectively: 


c ,  v/  r 

CseM  =  — 

n 


42 


(E  -  Er)1/2  H(E  -  Er 


(2.3) 


and 


„  42 

Csh(E)  =  ~ 
n 


i  *V/2 


mh 


\^2  j 


(Ev  -  E)1/2  H(EV  -  E) 


(2.4) 


where  H  (x)  is  the  Heavy-side  function  defined  as: 

fl  x  >  0 

H(x)  =  l 

0  otherwise 


(2.5) 


and  Ec,  Ev  are  conduction  and  valance  band  edges  defined  on 
Figure  2-2,  m*e ,  m*h  are  the  effective  masses  for  electrons 

and  holes  at  those  band  edges,  and  fi  is  short-hand  notation 
Planck's  constant  divided  by  2 n  . 

The  density  of  electron  concentration  in  electron 
energy  space  is  then  obtained  from  combining  Eqs .  (2.1)  and 

(2.2)  : 

dn 

Ce(E)  =  —2.  =  £se(E)  *  FD(E)  (2.6) 

dE 

and,  similarly,  for  holes  using  (2.2)  and  (2.3) : 

4(E)  =  =  4h(E)  *  FD(E)  (2.7) 

dE 

which  is  the  density  of  hole-concentration  in  electron 
energy  space. 

The  equilibrium  concentrations  of  electrons  and  holes 
are  compactly  expressed  as: 
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n0(EF)  =  J  Ce(E)  dE 


and 


p0(EF)  =  |  Ch(E)dE  (2.9) 


where  implicit  dependence  on  Fermi  energy  level  is  apparent 
through  substitutions  of  Eqs .  (2.1)  through  (2.7),  into 

Eq.  (2.8)  .  The  transformations,  set  as  x  =  (E-Ec)/kT  for 
evaluation  of  n0(EF)  and  set  as  x  =  (Ev-E)/kT  for  evaluation 
of  p0  (Ef)  ,  leads  to: 


na(EF) 


and 


where 


4~2 

r 

m 

*kT 'l 

e 

3/2 

v/r 

"  2  ° 

1/2  7 
x  dx 

7T2 

V 

h2  J 

2 

1  +  e  ' 

V2 

f 

m 

3/2 

2 

■  x1/z  dx 

n 2 

V 

h2  V 

2 

Jrt 

1  +  ex_r 

V 

=  (Ef 

_ 

Ec)/ 

kT 

(2 .10) 


(2.11) 


(2 . 12) 


and 

r  =  (Ev  -  Ef)  /  kT  .  (2.13) 

The  expressions  defined  by  (2.10)  and  (2.11)  can  be 
simplified  using  a  standard  approximation  [6]  as  shown  in 
the  next  paragraph. 


The  Fermi  integral  is  given  by  [6] : 


E1/2(y) 


2  r  x1 7  2dx 
~4n  ]  1  +  e<x_y> 


(2.14) 
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For  y  <  - 3: 


f 


f1/2(y) 


2 


J  e  Xx1 '  “dx 

V  v/‘  c  j 


-  e 1 


(2 . 15) 


This  approximation  applied  in  the  context  of  evaluation  of 
carrier  concentrations,  is  known  as  the  Boltzman' s  Approx¬ 
imation.  In  physical  terms,  the  applicability  of  the 
approximation  depends  on  the  Fermi  level  being  located  in 
the  bandgap  and  not  within  3kT  of  the  conduction  and 
valance  band  edges. 

After  combining  (2.10)  thru  (2.15)  it  follows  that: 

n0(EF)  =  Nce(EF~Ec)/ kT  Xn  (2.16) 


and 


p0(EF)  =  Nve^-E*)/kTXp 


(2.17) 


where  effective  electron  and  hole  density  of  states  are 
respectively : 

/  *  \  3  /  2 

(  2 nm„  kT) 


Nc  =  2 


and 


Nv  =  2 


V  ^  7 


(  27tm*kT^ 

v  *2  J 


(2 .18) 


3/2 


(2 .19) 


and  Xn  and  are  degeneracy  factors  described  in  the  next 
paragraph . 

The  electron  and  hole  concentration  degeneracy  factors 
[3]  are: 

%n  =  F1/2(ri)  /  e"  (2.20) 


for  which  X„~^  for  rl<~ 3  and: 

=  F1/2(y)  /  er 
11 


(2.21) 


for  which  ;jf  =lfor  y<- 3. 

If  the  Boltzman' s  Approximation  can  be  applied,  the 
usual  case  if  the  doping  is  not  too  heavy,  then  %n~ 1  and/or 
XP~^-  The  standard  Boltzman' s  approximation  form  for  the 

equilibrium  carrier  concentrations  is  recovered.  Speci¬ 
fically 

n0  «  Nce{EF~Ec)/kT  (2.22) 


and 


Equations 

intrinsic 


Po  “ 

(2.22)  and  (2.23) 
conditions,  i.e.. 


Nve(Ev~EP/kT  .  (2.23) 

can  be  taken  as  exact  under 
no  impurity  dopant  employed. 


Two  rules  are  needed  to  define  the  intrinsic 
condition.  These  rules  predict  the  Fermi  level  under 
intrinsic  conditions,  EF±r  and  the  extrinsic  carrier 
concentration,  n±.  Under  intrinsic  conditions,  electron  and 
holes  are  created  in  pairs  and  therefore  ni  =  no  =  po.  The 
application  using  (2.22)  and  (2.23)  leads  to: 

n±  =  =  ylNcNve-(*-B’)/2kT  (2.24) 


and  after  setting  EF  to  EF±  and  solving,  one  determines: 

(■ Ec  +  Ev'>  , 


kT  in 


\NcJ 


+ 


2 


2.251 


Equations  (2.24)  and  (2.25)  can  be  expressed  in  term 
of  the  bandgap,  Eg,  which  is  represented  on  Figure  2-2  as 
the  energy  difference  between  the  conduction  and  valence 
band.  If  the  valance  band  is  used  as  a  reference  energy 
level,  i.e.,  Er  = 0  as  represented  on  Figure  2-2,  Eg  is  then: 
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(2.26) 


_  Ec  Ev 
E’ - - 


Ec  +  Ev 
2 


and  it  follows  [6]  that  Eqs .  (2.22)  and  (2.23)  can  be 

expressed  in  terms  of  n±  and  EF±  as: 


Po 


n.e 


(Efi-Ef)/  kT 


(2.27) 


and 


n 


—  n.e 


(ef 


i)  / kT 


(2.28) 


Equations.  (2.27)  and  (2.28)  are  the  equilibrium  carrier 
concen-trat ions . 


Substitutional  impurities  can  exist  in  one  of  two 
states.  Specifically,  either  the  neutral  state  or  an  ion¬ 
ized  state  can  exist.  The  basic  rules  that  govern  donor 
and  acceptor  ionizations  [6]  are: 


IV, 


N, 


(2.29) 


and 


(2.30) 


where  subscript  ' d'  is  for  donor  and  ' a '  is  for  acceptor, 

Nd°  and  Na°  refer  to  neutral  donor  and  acceptor  impurity 
concentrations  respectively,  Nd+  and  Na~  refer  to  ionized 
donor  and  acceptor  concentrations,  and  gd  and  ga  are  occu¬ 
pation  indices.  There  are  electron  spin-based  arguments 
[6,12]  that  indicate  in  most  semiconductors  ga=  4  and  gd  =  2. 
The  fraction  of  ionization  can  be  established  via  several 
observations.  First,  the  sum  of  ionized  parts,  plus  neu¬ 
tral  must  be  added  to  the  total  impurity  used,  specifically 
for  donors : 
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(2.31) 


Nd  = 


Nd  +  N 


d 


and  for  acceptors: 

Na  =  N°a  +  N~  . 


(2.32) 


After  dividing  both  sides  of  Eq.  (2.31)  by  Ndr  it  is 
found  that : 


f 

V 


N 


d  J 


1  + 


K 

N* 


(2 . 33) 


where  Nd  and  Na  are  identified  on  Figure  2-2. 


Similarly,  after  dividing  both  sides  of  Eq.  (2.32)  by 
Na,  it  is  found  that: 


f 


V 


K 

N 


a  J 


l  + 


K 

n: 


(2.34) 


Therefore,  substitution  of  Eqs .  (2.31  and  2.32)  into  Eqs . 

(2.33  and  2.34)  leads  to: 


N 

N 


d_ 

d 


1 


1  +  9de 


(2 . 35) 


and 


1 


1  +  9ae 


(2 .36) 


which  provides  the  guidelines  for  Fermi  energy  level 
dependent  degrees  of  ionization. 

Donor  energy  levels  are  typically  cited  using  the 
conduction  band  as  the  reference.  For  similar  reasons,  the 
acceptor  energy  levels  are  typically  cited  using  the 
valence  band  as  a  reference  [3] : 

Ed  =  Ec  -  A Ed  (2.37) 


14 


and 


Ea  =  A  Ea  (2.38) 

which  is  consistent  with  Figure  2-2. 

Assuming  we  have  Md,  mono-valent  donors,  and  Ma,  mono¬ 
valent  acceptors,  the  conditions  of  local  charge  neutrality 
exist.  This  is  based  on  assumptions  #1  and  #2.  For  the 
condition : 

Md  Ma 

P{T,Ef)  =  Po  -  n0  +  ~  SK"),  =  0  (2.39) 

1=1  1=1 

where  there  is  an  implicit  dependence  on  both  EF  and  the 
temperature  within  each  term  of  Eq.  (2-22) .  The 
dependencies  are  obtained  from  Eqs .  (2.35  to  2.38)  .  For  a 
specified  temperature,  Eq.  (2.39)  results  in  a  non-linear 

expression  in  Z,  where  Z  =  e  P1  F  At  a  specified  Z,  the 

carrier  concentrations  can  be  predicted  from: 


and 

po  =  n±Z  (2.41) 

which  is  an  equivalent  revised  version  of  Eqs.  (2.27)  and 
(2.28)  . 

Making  use  of  the  Fermi  probability  distributions  to 
predict  ionization  levels  of  the  donors  and  acceptors, 
specifically  this  requires  the  substitution  of  Eqs.  (2.35, 
2.36,  2.40,  and  2.41)  into  Eq.  (2.39)  and  dividing  by  the 
intrinsic  concentration  which  leads  to  a  condition  on 
P  (Z,T)  : 

A  Ndl  A  Nai 

P(Z)  =  o  =  Z  +  y  - - — —  -  Z  -  y  - - -  .  (2.42) 

ti  1  +  V1  ti  1  +  KaZ 
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The  over-bar  notation  indicates  division  by  n±.  The  tempera¬ 
ture  dependent  ionization  factors: 


K 


(2 .43) 


and 


(2.44) 


serves  as  a  gauge  for  degrees  of  ionization  of  the  associ¬ 
ated  impurity.  In  summary,  Eqs .  (2.42  to  2.44)  serve  as  a 

mathematical  basis  for  both  the  analytic  and  numerical 
[3,4]  MATLAB-based  approaches. 


C.  ANALYTIC  SOLUTION  APPROACH 

This  situation  is  described  by  the  special  case  that 
Ma  = 1  and  Md  = 0  .  The  summation  subscript  in  Eq.  (2.42)  can 
be  dropped  to  lighten  the  notation.  Straightforward  algebra 
leads  to  a  cubic  equation  in  Z : 


Z3  +  Z2 
if. 


Na 


V1  +  K*J 


—  =  0  . 

if. 


( 2 .45] 


The  corresponding  result  for  a  single  impurity  N-type 
is  exactly  the  same  form,  producing  a  cubic  in  Y-\!Z  with 
revised  coefficients  obtained  by  letting  the  ' a '  subscript 
(for  acceptor)  be  replaced  with  the  '  d'  subscript  (for 
donor) .  The  solution  can  be  defined  in  terms  of 
coefficients  for  the  reference  cubic  equation: 

x3  +  azx2  +  azx  +  a3  =  0  .  (2.46) 


To  facilitate  representation  of  a  solution,  the 
following  intermediate  parameterization  of  the  problem  is 
commonly  taken  [9] : 

Q  =  (3a2  -  a3)  /  9  (2.47) 
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and 


R  =  (93,32  -  27a3  -  2a,3)/  54  (2.48) 

and 


D  =  Q3  +  R2 


(2.49) 


where,  D,  which  is  referred  to  as  the  "discriminant"  for 
the  cubic  problem,  will  dictate  the  type  of  solutions 
possible.  As  per  the  fundamental  theorem  of  algebra  there 
will  be  three  roots.  For  D  negative,  all  roots  are  real 
while  for  D  positive  only  one  root  is  real  while  the  other 
two  are  complex  conjugates.  If  D  is  zero  there  will  be 
repeated  root.  For  the  problem  being  evaluated  it  can  be 
shown  using  symbolic  mathematical  methods  that  Q  and  D  will 
always  be  negative.  Furthermore,  it  can  then  be  shown  that 
of  the  three  roots,  only  one  is  positive  and,  therefore, 
physically  acceptable.  That  root  will  be  predicted  by  the 
following  formula  [9] : 


Z 


x 


2 4~Q  cos  (i 9  /  3) 


(2.50) 


where  0  -  cos  1  i J-Q 3  j  once  Z  is  determined,  the  majority 
carrier  concentration,  p0,  can  be  predicted  from  Eg.  (2.41)  . 


D .  CHAPTER  SUMMARY 

This  Chapter  has  established  the  mathematical  basis 
for  both  MATLAB-based  analytical  and  numerical  predictor 
algorithms  for  temperature  dependence  of  the  majority 
carrier  (see  Figure  2-1) .  In  this  Chapter  the  main  steps 
leading  to  a  single  equation  analytic  solution  are 
described.  In  the  next  Chapter,  the  third  approach,  based 
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on  the  commercially  available  Silvaco  International  Soft¬ 
ware,  is  briefly  explained. 
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III.  SIMULATION 


A.  CHAPTER  INTRODUCTION 

For  this  project,  a  suite  of  commercially  available 
simulation  tools  [1]  intended  for  semiconductor  testing  was 
used.  The  primary  objective  of  the  thesis  was  to  generate 
Silvaco  predictions  for  the  temperature  dependent  majority 
carrier.  A  comparison  will  be  made  with  the  MATLAB-based 
analytical  and  numerical  methods  discussed  in  the  previous 
chapter  (see  Figure  2-1) .  The  data  from  all  three 
approaches  were  then  compared  graphically  using  Microsoft 
Excel.  The  following  sections  will  discuss  the  processes 
developed  in  the  simulation  to  generate  the  data  for  this 
comparison  testing. 


B.  SILVACO  SIMULATION  TOOLS 

The  simulation  in  Silvaco  was  designed  using  the 
following  modules  from  the  suite:  ATLAS,  DeckBuild, 
TonyPlot  and  DevEdit .  The  following  outline  (Figure  3-1) 
illustrates  the  interconnecting  relationships  between  the 
individual  modules. 
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Figure  3-1.  Silvaco  Interaction  Flowchart  (From 

[1]  )  • 

ATLAS  is  a  2D/3D  simulator  for  semiconductor  devices. 
ATLAS  can  provide  data  and  insights  into  the  internal 
physical  mechanisms  of  a  device  based  on  predicted 
electrical  behavior.  It  can  either  be  used  as  a  stand¬ 
alone  tool  or  as  a  core  unit  for  the  Silvaco  Virtual  Wafer 
Fabrication  environment.  This  module  was  used  in 
connection  with  the  other  modules  to  deliver  powerful  and 
accurate  data  on  the  behavior  of  semiconductors  [1] . 

DeckBuild  is  an  interactive  graphical  user  interface 
to  provide  a  user-friendly  runtime  environment  for 
integration  of  the  different  aspects  of  the  Silvaco 
Software  suite.  A  control  window  is  provided  for  file 
creation  and  control.  Many  of  the  features  are  automated 
to  allow  for  accurate  simulation  in  a  simple  to  use 
environment.  TonyPlot  is  the  stand-alone  program,  which 
can  also  be  referenced  in  DeckBuild,  to  display  the  results 
[1]  • 

DevEdit  is  a  device  editor  which  can  be  used  to 
generate  a  mesh  for  the  structure  designed  in  DeckBuild  and 
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ATLAS.  A  limitation  of  device  simulators  prior  to  DevEdit 
was  inadequate  or  poor  structure  meshes.  DevEdit' s  usage 
was  integrated  into  DeckBuild  to  allow  for  a  more  complete 
and  accurate  solution  [1] . 

C.  DEVELOPMENT  OF  THE  SILVACO  SIMULATION  PROCESS 

In  order  to  develop  the  overall  simulation,  the  first 
step  was  to  understand  the  full  scope  of  the  problem  and  to 
determine  the  best  features  of  Silvaco  to  find  a  solution. 
There  were  many  automated  features  and  defaults  settings 
needed  to  be  addressed  to  ensure  that  both  sets  of  data 
were  based  on  the  same  fundamentals.  The  DeckBuild  coding 
in  Silvaco  followed  the  standard  format  of  the  program  as 
illustrated  in  Figure  3-2: 


Group 

Statements 

1.  Structure  Specification 

MESH 

REGION 

ELECTRODE 

DOPING 

2.  Material  Models  Specification 

MATERIAL 

MODELS 

CONTACTS 

INTERFACE 

3.  Numerical  Models  Specification  - 

METHOD 

4.  SnlnHnn  Sppr.i  fi  rail  nn 

LOG 

SOLVE 

LOAD 

SAVE 

5.  Results  Analysis 

EXTRACT 

TONYPLOT 

Figure  3-2.  Input  File  Layout  (After  [1]) . 
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In  Figure  3-2,  the  left  side,  listed  under  Group,  outlines 
how  commands  are  clustered  in  Silvaco.  The  right  side, 
listed  under  Statements,  catalogs  several  of  the  commonly 
used  commands  and  the  order  in  which  they  are  called. 


Several  important  values  were  set  at  the  beginning  of 
the  DeckBuild  code  to  ensure  that  default  settings  were  not 
used.  The  specified  parameters  included:  the  dopant  used, 
the  e-alpha  (a)  ,  e-beta  ((3)  ,  and  e-band  gap  {Ego)  values, 
which  were  set  before  initializing  the  program.  As  pre¬ 
viously  stated,  the  bandgap  represents  one  of  the  most 
important  parameters  of  a  semiconductor  and  thus,  it  was 
set  explicitly  in  DeckBuild  coding.  The  equation  for 
bandgap  narrowing  listed  in  the  Silvaco  ATLAS  Manual  [1]  is 
as  follows: 


E  (T) 

g 


Eg  (3  0  0)  +  a 


30  02 

300  +  /3 


rj-,2 

T  +  j8 


(3.1) 


where  T  is  the  temperature  in  degrees  Kelvin  and  the  a  and  (3 
parameters  are  semiconductor  material  dependent.  Eg( 300)  is 
the  energy  gap  at  300°K.  An  equivalent  f ormulization  used 
in  the  MATLAB-based  program  is  in  terms  of  Eg(0)  instead  of 
£g  (300)  [6,  10]. 

The  structural  dimensions  of  the  device  were  then  set 
and  the  contacts  were  placed  with  DevEdit .  The  device 
material  was  then  set  and  the  doping  concentration  was  also 
set  in  DeckBuild.  Choosing  a  model  was  the  next  challeng¬ 
ing  aspect  of  research  since  certain  DeckBuild  coding 
options  would  decide  the  included  physical  effects.  This 
was  needed  in  order  to  obtain  the  most  comprehensive 
solution  to  compare  with  the  MATLAB-based  approaches  (see 
Figure  2-1 ) . 
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For  the  model  options  in  the  Silvaco  Simulations 
DeckBuild  code,  incomplete,  conmob,  and  fermi  were  chosen, 
because  they  covered  the  widest  comprehensive  range  of 
temperature-dependent  effects.  The  incomplete  portion  was 
used  to  properly  account  for  freeze-out  when  there  were  not 
sufficient  amounts  of  thermal  energy  to  assume  100%  dopant 
ionization.  The  conmob  function  was  used  to  maintain 
concentration  mobility  as  per  the  simple  power  rule  for 
temperature  dependence.  The  fermi  command  was  used  to 
account  for  cases  of  heavily  doped  materials. 

Along  with  the  technical  challenges  of  understanding 
the  correct  configuration  for  modeling  a  semiconductor, 
there  were  the  difficulties  in  manipulating  Silvaco  due  to 
inherent  limitations  of  the  program  suite.  A  major 
obstacle  related  to  the  inability  to  program  a  temperature 
sweep  in  Silvaco  within  a  single  run.  A  method  had  to  be 
developed  to  complete  this  process  in  an  automated  fashion. 
It  was  accomplished  by  replicating  single-run  code  segments 
within  DeckBuild.  Each  of  the  code  segments  was  essen¬ 
tially  the  same,  except  for  the  temperature,  which  was 
incremented  consecutively  (see  Appendix  B) .  The  tempera¬ 
ture  along  with  other  pertinent  data  were  recorded  and 
logged  after  each  step,  which  completed  the  ATLAS  run. 

This  process  was  repeated  over  a  selected  temperature 
range.  Stored  data  could  then  be  exported  to  a  program, 
such  as  Microsoft  Excel.  This  process  allowed  for 
significant  flexibility  in  selecting  the  semiconductor 
parameters.  It  also  allowed  for  an  automation  in  the 
process  which  required  very  minimal  user  input.  The 
overall  process  is  illustrated  in  Figure  3-3. 
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Figure  3-3.  Flowchart  of  Silvaco  Process. 

D.  DEVELOPMENT  OF  THE  MOVIE 

TonyPlot  was  used  to  create  a  slideshow,  to  aid  in 
visualization  of  the  process.  The  structure  files  from  the 
DeckBuild  runs  were  saved  and  TonyPlot  was  used  to  display 
the  hole/electron  concentrations  of  the  lattice  over  the 
temperature  steps.  These  data  outputs  were  then  exported 
into  other  commercial  software  packages  that  allowed  for 
graphic  file  animation,  such  as  Microsoft  Office  Power¬ 
Point.  Two  samples  of  the  movie  are  provided  in  Figure  3-4 
and  3-5.  The  arrow  over  the  curves  on  the  right  side  iden¬ 
tifies  the  temperature.  The  upper  half  shows  a  selected 


movie  slide  when  the  temperature  was  approximately  10  K 
while  the  lower  half  was  taken  at  a  much  higher  temperature 
of  900  K.  Note  that  colors  for  the  left  side  are  legend- 
encoded  for  the  concentrations  of  the  majority  carrier  at 
the  corresponding  temperature.  More  details  on  the 
procedure  for  creating  a  movie  are  provided  in  Appendix  B. 


E .  CHAPTER  SUMMARY 

This  Chapter  outlined  and  explained  the  technical 
aspects  of  utilizing  commercially  available  Silvaco 
International  software  for  the  objectives  of  this  thesis. 

A  more  detailed  description  of  this  procedure  is  provided 
in  Appendix  B.  The  next  chapter  will  provide  a  review  and 
comparison  of  the  Silvaco  based  results  with  predictions 
generated  from  MATLAB-based  analytic  and  numerical 
algorithm. 
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Plot  of  3-regimes  of  Silicon 


TenpeutiieiK) 


Figure  3-4 : 


Selected  Samples  of  50  slide  Movie  T  = 


60°K 
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Plot  of  3-regimes  of  Silicon 


TenpeutiieiK) 


Figure  3-5: 


Selected  Samples  of  50  slide  Movie  T  = 
1 0  0  0°K 


27 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


28 


IV.  RESULTS  ANALYSIS  AND  COMPARISON 


A.  CHAPTER  INTRODUCTION 

This  chapter  provides  a  comparison  of  the  results  for 
various  approaches  obtained  from  examples  taken  of  Silicon 
and  GaAs .  The  different  conditions  and  cases  documented 
will  be  discussed.  For  this  project,  Silvaco  Suite  was 
used  to  develop  a  procedure  to  predict  temperature 
dependence  of  majority  carrier  concentration.  This  was 
discussed  in  some  detail  in  the  previous  chapter.  The 
results  from  Silvaco  are  compared  with  MATLAB-based 
analytic  and  numerical  predictions.  The  data  from  both 
simulations  are  then  compared  graphically  in  Microsoft 
Excel.  The  following  sections  will  discuss  the  processes 
developed  in  the  simulation  to  generate  the  data  for  this 
comparison  testing. 


B.  RESULTS  FOR  SILICON 


Figure  4-1  contains  six  plots  of  logio  (Po)  versus 
1000* (1/T)  for  various  conditions  from  the  various  sources 
(numerical,  analytic,  and  Silvaco) .  The  data  was  also 
grouped  into  the  case  with  and  without  bandgap  narrowing 
(BGN) .  The  model  for  bandgap  narrowing  [6]  used  here  is: 


(XT 

E  (T)  =  E - 

g  g  (T  +  p) 


(4.1) 


Numbers  codes  for  the  various  curves  on  Figure  4-1  appear 
with  the  description  in  the  legend.  Recommended  values 
[2,15]  for  silicon,  including  BGN  parameters  a  and  (3,  appear 
in  Table  4-1.  Note,  as  expected,  curves  from  the  MATLAB- 
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based  numerical  and  analytic  methods  are  indistinguishable. 
For  example,  MATLAB-based  analytic  curve  #4 (w/BGN)  overlays 
MATLAB-based  numerical  curve  #6 (w/BGN).  Furthermore,  the 
impact  of  BGN  is  apparent  in  the  curves  at  the  higher 
temperatures . 


With  and  Without  Bandgap  Narrowing  Plot 


Figure  4-1  .  Composite  Log  Plot  of  three  Regimes  of 

Silicon . 
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The  difference  between  Silvaco  predictions  and  the 
one-dimensional  modeling  efforts  are  slight  but  noted  for 
future  investigation. 


Parameter 

Value 

Ea 

0.045  eV 

Na 

1017  I/cm3 

Ego 

1.11  eV  (silicon) 

9a 

4 

BGN  a 

0 . 00047 

BGN  /3 

636.0 

Table  4-1.  Silicon  Parameters. 


C.  RESULTS  FOR  GALLIUM  ARSENIDE 

The  results  for  GaAs  seen  in  Figure  4-2  were  similar 
to  results  appearing  in  Figure  4-1.  The  parameters  used 
for  this  test  are  provided  in  Table  4-2  [2,  15] : 


Parameter 

Value 

Ea 

0.15  eV 

Na 

1012  I/cm3 

Ego 

1.42  eV  (GaAs) 

9a 

4 

BGN  a 

0 . 00054 

BGN  /3 

204 . 0 

Table  4-2.  GaAs  Parameters. 
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With  and  Without  Bandgap  Narrowing  Plot  for  GaAs 


1000  X  1/T(K) 


Figure  4-2.  Composite  Log  Plot  of  three  Regimes  of 

GaAs  . 

Again,  all  three  regimes,  i.e.,  freeze-out,  exhaustion,  and 
intrinsic,  were  evident  on  Figure  4-2.  As  before,  in  the 
exhaustion  regime,  the  majority  carrier  was  approximately 
equal  to  the  doping  level,  i.e.,  1012  #/cm3.  The  curves 
indicated  that  the  MATLAB-based  approaches  produced  pre¬ 
dictions  that  were  indistinguishable.  It  should  also  be 
noted  that  the  level  of  agreement  between  the  Silvaco-based 
prediction  and  the  MATLAB-based  approaches  was  better  in 

the  GaAs  case  (see  Figure  4-2)  versus  the  silicon  case  (see 
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Figure  4-1) .  This  is  especially  true  in  comparing  the 
level  of  agreement  in  the  freeze-out  regimes  of  both 
Figures  4-1  and  4-2.  As  it  turns  out,  the  testing  pro¬ 
cedure  for  the  two  materials.  Si  and  GaAs,  was  essentially 
the  same  with  one  important  exception.  For  the  Silvaco- 
based  simulation  of  Silicon,  the  input  specification  for 
the  acceptor  energy  level  (Ea)  was  not  done  directly. 

Instead,  the  chemical  impurity  dopant  Boron  was  specified 
in  the  DeckBuild  code.  In  the  case  of  the  GaAs  test,  this 
feature  was  overridden  by  a  direct  specification  of  the 
acceptor  energy  level  ( Ea )  in  the  DeckBuild  code.  The 
latter  approach  is  clearly  aligned  with  the  MATLAB-based 
methods.  It  is  a  reasonable  conjecture  that  this  accounts 
for  the  noticeable  improvement  in  the  level  of  agreement 
between  the  Silvaco  and  MATLAB-based  methods. 

D.  DISCUSSION  OF  RESULTS 

The  level  of  agreement  between  results  generated  from 
the  MATLAB  approaches  and  from  the  Silvaco  source  code 
tended  to  confirm  the  accuracy  of  the  models.  The  Silvaco 
and  MATLAB-based  approaches  were  in  close  agreement  in 
predictions  for  the  majority  the  carrier  concentration  over 
all  three  single  impurity  regimes  of  freeze-out,  exhaus¬ 
tion,  and  intrinsic. 

An  important  effect  included  in  the  testing  was 
Bandgap  Narrowing  (BGN)  [6]  .  The  results  were  categorized 
for  each  material  tested  with  BGN  and  without  BGN.  The 
differences  observed  for  BGN  proved  to  be  fairly 
significant  at  higher  temperatures.  In  the  exhaustion  and 
freeze-out  regimes,  the  impact  of  BGN  was  not  noticeable. 
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E .  CHAPTER  SUMMARY 

Both  the  one-dimensional  analytic  model  and  numerical 
algorithm  were  essentially  indistinguishable  in  terms  of 
predictions.  As  expected,  when  BGN  effects  were  included 
in  the  models,  the  impact  was  observed  at  the  higher 
temperatures.  The  difference  between  predictions  of  the 
Silvaco-based  two-dimensional  modeling  and  the  one¬ 
dimensional  solutions  were  observed  to  be  slight  but  worthy 
of  future  investigation.  The  next  Chapter  will  summarize 
the  main  conclusions  and  will  provide  suggestions  for 
future  work. 
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V. 


CONCLUSIONS  AND  RECOMMENDATIONS 


A.  SUMMARY  OF  WORK 

This  Chapter  provides  a  summary  of  the  thesis  and 
provides  some  final  thoughts  on  future  areas  of  study  and 
improvement.  The  principal  completed  goal  of  this  research 
has  been  to  design  a  convenient  automated  procedure  that 
will  predict,  via  Silvaco  International  semiconductor 
simulation  software  [1],  temperature-dependent  exhaustion 
and  freeze-out  effects  in  single-dopant,  homojunction 
semiconductors .  To  gauge  the  level  of  success  of  this 
effort  the  results  generated  from  the  Silvaco  simulation 
software  have  been  compared  with  similar  predictions  from 
MATLAB-based  analytical  and  numerical  algorithms.  This  was 
done  for  both  silicon  and  gallium  arsenide  semiconductors. 

Chapter  I  provided  an  introduction  to  the  objectives 
and  benefits  of  this  research,  along  with  a  preview  of  what 
topics  would  be  addressed.  Chapter  II  presented  an  in- 
depth  mathematical  analysis,  which  serves  as  a  basis  for 
MATLAB-based  analytic  and  numerical  approach.  This  chapter 
also  detailed  the  process  used  to  determine  the  analytical 
solution.  Chapter  III  was  focused  on  imparting  an  overview 
of  Silvaco  and  how  it  was  implemented  to  provide  the  neces¬ 
sary  results.  Chapter  IV  was  used  to  present  the  data 
acquired  from  the  different  tests  and  to  make  a  comparative 
study  based  on  the  data.  A  sample  selection  of  slides  from 
the  Silvaco  movie  prepared  during  the  thesis  is  also 
presented  in  Chapter  IV. 
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B.  CONCLUSION 

This  thesis  was  designed  to  be  part  of  a  larger  study 
in  creating  very  robust  MATLAB-based  analytic  and  numerical 
simulations  to  provide  an  accurate  portrait  of  temperature 
effects  occurring  in  a  semiconductor.  This  would  allow  for 
a  versatile  and  reliable  method  in  predicting  the  regimes 
of  a  semiconductor. 

After  reviewing  the  data,  it  can  be  concluded  that  the 
MATLAB-based  analytical  and  numerical  models  compared  very 
favorable  to  the  results  obtained  in  Silvaco.  All  three 
methods  accurately  predict  the  semiconductor  behavior  in 
the  three  regions  from  freeze-out  to  the  intrinsic.  The 
MATLAB-based  analytic  and  numerical  coding  allowed  for  a 
quick  and  robust  method  to  generate  the  majority  carrier 
concentration  data  of  a  semiconductor  under  a  varying 
degree  of  changes  in  material  and  doping  level.  Currently, 
the  computation  cost  of  results  obtained  via  MATLAB 
required  much  less  time  than  the  Silvaco  method  by  roughly 
a  factor  of  50. 

C .  RECOMMENDATIONS 

There  has  been  significant  progress  made  in  this 
research,  but  there  still  exist  much  additional  work  to  be 
accomplished  in  this  area.  The  future  of  analytical 
modeling  in  MATLAB  needs  to  be  refined  and  expanded  to 
allow  for  a  greater  range  of  effects  and  conditions.  For 
example,  the  range  of  cases  that  can  be  solved  analytically 
could  be  widened  to  include  multiple  impurities.  This 
capability  already  exists  in  the  MATLAB-based  numerical 
algorithm.  At  this  present  time,  the  numerical  algorithm 
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has  not  satisfactorily  been  tested  for  high  levels  of 
doping  when  the  Boltzman' s  approximation  fails. 

There  still  exist  many  conditions  and  parameters  in 
Silvaco  which  have  not  been  accounted  for  which  could  be 
implemented  in  a  faster  MATLAB  model.  A  more  efficient 
method  of  handling  temperature  ranges  in  Silvaco  is  also 
something  that  needs  to  be  improved  upon. 

Finally,  utilizing  Silvaco  to  create  and  test  a  three- 
dimensional  model  might  provide  more  conclusive  data  to 
continue  to  validate  the  results  obtained  for  the 
analytical  model.  This  would  provide  a  significant  insight 
into  the  limits  associated  with  the  one-dimensional  models 
being  tested  in  MATLAB. 
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APPENDIX  A.  INTRODUCTION  TO  GENERAL  SEMICONDUCTOR 
PRINCIPLES  (AFTER  [11,  12]) 


A.  QUANTUM  MECHANICS 

To  fully  understand  the  basis  of  this  thesis,  a  basic 
knowledge  of  quantum  effects  is  required  to  understand  why 
particles  behave  as  they  do  at  the  sub-atomic  and  atomic 
dimensional  scale.  Before  the  twentieth  century,  scien¬ 
tists  had  a  firm  belief  that  the  physical  world  was  divided 
into  two  worlds,  the  world  of  particles  and  the  world  of 
waves.  This  belief  was  based  on  the  classical  system 
theorized  by  Sir  Isaac  Newton. 

In  the  20th  Century,  this  particle/wave  paradigm  was 
challenged  and  proven  to  be  incomplete,  based  on  the  early 
work  of  physicist  such  as  Neil  Bohr.  Werner  Heisenberg, 
and  an  Austrian,  Erwin  Schrodinger  postulated  a  new  Quantum 
Theory  based  on  the  aforementioned  work,  and  thus  a  valid 
and  proven  method  was  constructed  to  provide  a  means  of 
describing  sub-atomic  particles. 

Quantum  mechanics  theorizes  that  the  electronic  struc¬ 
ture  of  an  atom  is  based  on  quantum  conditions.  These 
conditions  are  postulated  in  the  theory  that  electrons  in 
an  atom  are  restricted  to  certain  energy  level.  This 
"quantization"  defines  certain  allowable  transitions  invol¬ 
ving  absorption  and  emission  of  energy  by  the  electrons 
[11]  • 

Neils  Bohr  determined  that  electrons  in  orbits  moved 
with  well  defined  momentum.  Based  on  this  work  and  through 
absorption  and  emission,  an  electron  may  transition  to  the 
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next  smaller  or  larger  energy  level  based  on  the  following 
rules  : 

1)  Electrons  exist  in  certain  stable,  circular 
orbits  about  the  nucleus. 

2)  The  electron  may  shift  to  an  orbit  of  higher 
or  lower  energy,  thereby  losing  or  gaining 
energy  equal  to  the  difference  in  the  energy 
levels  (by  absorption  or  emission  of  a  photon 
of  energy  hv,  see  Figure  A-l)  . 

3)  The  angular  momentum  Pg  of  the  electron  in  an 
orbit  is  always  an  integral  multiple  of 
Planck's  constant  divided  by  2 n  [11]  . 


E2 

hv 


El 


Figure  A-l.  Emission  of  Photon. 


B.  SEMICONDUCTOR  MATERIAL  PROPERTIES 

A  semiconductor  is  able  to  conduct  electricity  at  room 
temperature  more  readily  than  an  insulator,  but  less  easily 
than  a  metal.  Electrical  conductivity,  which  is  the 
ability  to  conduct  electrical  current  under  the  application 
of  a  voltage,  has  one  of  the  widest  ranges  of  values  of  any 
physical  property  of  matter.  Such  metals  as  copper, 
silver,  and  aluminum  are  excellent  conductors,  but  such 
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insulators  as  diamond  and  glass  are  very  poor  conductors. 

At  low  temperatures,  pure  semiconductors  behave  like  insu¬ 
lators.  Under  higher  temperatures  or  light  or  with  the 
addition  of  impurities,  however,  the  conductivity  of 
semiconductors  can  be  increased  dramatically,  reaching 
levels  that  may  approach  those  of  metals  [12] . 

The  common  semiconductors  include  chemical  elements 
and  compounds  such  as  silicon,  germanium,  gallium  arsenide, 
and  lead  telluride.  They  increase  in  conductivity  with 
temperature,  light,  or  impurities  arising  from  an  increase 
in  the  number  of  conduction  electrons,  which  are  the 
carriers  of  the  electrical  current.  In  a  pure,  or  intrin¬ 
sic,  semiconductor  such  as  silicon,  the  valence  electrons, 
or  outer  electrons,  of  an  atom  are  paired  and  shared  bet¬ 
ween  atoms  to  make  a  covalent  bond  that  holds  the  crystal 
together.  These  valence  electrons  are  not  free  to  carry 
electrical  current.  To  produce  conduction  electrons, 
temperature  or  light  is  used  to  excite  the  valence  elec¬ 
trons  out  of  their  bonds,  leaving  them  free  to  conduct 
current.  Deficiencies,  or  "holes,"  are  left  behind  that 
contribute  to  the  flow  of  electricity.  (These  holes  are 
said  to  be  carriers  of  positive  electricity.)  This  is  the 
physical  origin  of  the  increase  in  the  electrical  conduc¬ 
tivity  of  semiconductors  with  temperature.  The  energy 
required  to  excite  the  electron  and  hole  is  called  the 
energy  gap. 

For  a  semiconductor,  this  energy  gap  is  termed  the 
bandgap .  It  is  a  forbidden  region  in  which  an  electron 
cannot  exist,  as  it  is  the  "space"  between  the  quantized 
states.  The  upper  region  is  called  the  conduction  band  and 
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the  lower  is  labeled  the  valence  band.  This  separation 
between  the  energy  of  the  two  bands  called  the  bandgap  is 
labeled  as  Eg.  It  is  the  most  important  parameter  when 
dealing  with  semiconductors  since  it  is  the  basis  for  any 
utilization  of  a  given  material.  These  bandgaps  define  the 
relative  ease  to  "free"  electron  through  excitation  from 
thermal,  optical,  or  electrical  energy.  Semiconductors 
have  lower  bandgaps  (typically  0.5-4  eV)  while  insulators 
have  higher  bandgaps  (5-10  eV)  of  transition. 

Semiconductor  devices  have  many  varied  applications  in 
electrical  engineering.  Recent  engineering  developments 
have  yielded  small  semiconductor  chips  containing  millions 
and  in  the  near  future,  hundreds  of  millions,  of  transis¬ 
tors.  These  chips  have  made  possible  great  miniaturization 
of  electronic  devices.  More  efficient  use  of  such  chips 
has  been  developed  through  what  is  called  complementary 
metal-oxide  semiconductor  circuitry,  or  CMOS,  which  con¬ 
sists  of  pairs  of  p-  and  n-channel  transistors  controlled 
by  a  single  circuit. 


C .  CONDUCTION 

Certain  critical  terms  and  concepts  must  be  introduced 
to  better  understand  the  fundamental  of  semiconductors. 

Two  parameters  in  the  handling  and  usage  of  semiconductors 
are  temperature  and  doping.  These  two  factors  are  extreme¬ 
ly  important  in  determining  how  effect  a  semiconductor  is 
and  how  it  can  most  effectively  be  used  in  any  system 
design . 

Doping  is  a  method  to  produce  free  carriers  of  elec¬ 
tricity  in  a  semiconductor  by  adding  impurities  to,  or  to 
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"dope,"  the  semiconductor.  The  difference  in  the  number  of 
valence  electrons  between  the  doping  material,  or  dopant 
(either  donors  or  acceptors  of  electrons) ,  and  host  gives 
rise  to  negative  (n-type)  or  positive  (p-type)  carriers  of 
electricity.  This  concept  is  illustrated  in  the  accom¬ 
panying  diagram  of  a  doped  silicon  (Si)  crystal.  Each 
silicon  atom  has  four  valence  electrons  (represented  by 
dots  in  Figure  A-2);  two  are  required  to  form  a  covalent 
bond.  In  n-type  silicon,  atoms  such  as  phosphorus  (P)  with 
five  valence  electrons  replace  some  silicon  and  provide 
extra  negative  electrons.  In  p-type  silicon,  atoms  with 
three  valence  electrons  such  as  aluminum  (Al)  lead  to  a 
deficiency  of  electrons,  or  to  holes,  which  act  as  positive 
electrons.  The  extra  electrons  or  holes  can  conduct 
electricity . 


Figure  A-2 . 
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Si  Material  Doped  with  P  (From  [13]). 


When  p-type  and  n-type  semiconductor  regions  are 
adjacent  to  each  other,  they  form  a  semiconductor  diode, 
and  the  region  of  contact  is  called  a  p-n  junction.  (A 
diode  is  a  two-terminal  device  that  has  a  high  resistance 
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to  electric  current  in  one  direction  but  a  low  resistance 
in  the  other  direction)  The  conductance  properties  of  the 
p-n  junction  depend  on  the  direction  of  the  voltage,  which 
can,  in  turn,  be  used  to  control  the  electrical  nature  of 
the  device.  Series  of  such  junctions  are  used  to  make 
transistors  and  other  semiconductor  devices  such  as  solar 
cells,  p-n  junction  lasers,  rectifiers,  and  many  others. 

Temperature  is  also  very  important  and  as  it  is  raised 
from  0°K,  electrons  are  excited  from  the  thermal  energy  and 
they  move  across  the  bandgap  from  the  valence  band  to  the 
conduction  band.  The  resulting  material  is  one  with  some 
unoccupied  states  in  the  valence  band  and  this  vacancy 
creates  a  "hole",  or  an  empty  state.  This  creates  an 
electron-hole  pair  (EHP)  and  as  these  "freed"  electrons 
move  through  the  material,  this  motion  induces  an  electric 
field  and  current  flows. 

There  are  two  major  classes  of  semiconductors,  intrin¬ 
sic  and  extrinsic.  An  intrinsic  semiconductor  is  one  in 
which  the  electrical  properties  of  device  is  determined  by 
the  host  material,  such  as  with  bulk  silicon.  In  an 
extrinsic  semiconductor,  the  chemical  impurities  used  to 
dope  the  host  material  are  the  determining  factor  for  the 
electrical  properties  of  the  material.  These  dopants  can 
drastically  modify  the  electrical  properties  of  a  device  as 
with  gallium  arsenide. 
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As  illustrated  in  Figure  A-3,  when  the  temperature 
increases,  the  intrinsic  properties  of  the  material  over¬ 
ride  any  level  of  doping  applied  to  the  material.  This 
region  varies  based  on  material  and  doping  level  of  im¬ 
purities.  At  low  temperatures,  the  energy  from  thermal 
effects  is  not  high  enough  to  fully  activate  all  the  donor 
and  acceptor  impurities.  This  region  is  known  as  the 
freeze-out  region  and  this  topic  will  be  addressed  later  in 
this  report. 

If  temperature  increases  above  0  K,  excitation  of 
electrons  is  induced  and  EHPs  are  formed  as  thermally 
excited  electrons  cross  the  bandgap  to  the  conduction  band. 
The  positions  of  the  free  electron  and  hole  in  the  crystal 
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are  spread  out  over  several  lattice  spacing  and  can  be 
considered  quantum  mechanically  by  probability  distri¬ 
butions.  At  a  given  temperature,  there  is  a  certain  amount 
of  EHPs.  In  a  steady  state  condition  the  generation  rate 
(Gi)  of  EHPs  in  an  intrinsic  material  must  be  equal  to  the 
rate  of  recombination  of  electrons  and  holes  [11] . 

As  temperature  increases,  the  probability  distribution 
of  an  electron  in  a  range  of  allowed  energy  at  equilibrium 
follows  Fermi-Dirac  statistics,  because  they  obey  the  Pauli 
exclusion  principal.  The  Fermi-Dirac  distribution  function 
is  equal  to: 

/(£)  = - - -  (A— 1 ) 

where  Ep  is  called  the  Fermi  Level.  For  intrinsic  ma¬ 
terials,  the  probability  an  electron  is  at  the  Fermi  level 
is  very  near  U.  When  T=0K,  this  function  is  a  unit  step 
and,  as  T  increases,  the  function  takes  on  a  more  curved 
appearance  (see  Figure  A-4). 
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APPENDIX  B.  OPERATIONAL  GUIDE  TO  REPRODUCING  THE 

RESULTS  IN  SILVACO 


A.  OPERATIONAL  GUIDE  TO  REPRODUCING  SILVACO  DATA 

ACQUISITION  CAPABILITIES 

Silvaco' s  semiconductor  software  is  a  powerful  and 
very  capable  tool  in  deciphering  the  electrical  and 
material  properties  of  a  semiconductor.  The  vast  array  of 
built-in  functions  and  add-on  modules  create  a  powerful 
tool,  albeit  a  tool  with  a  steep  learning  curve.  This 
concise  section  will  attempt  to  help  create  a  standardize 
method  to  collect,  organize,  and  extract  the  data  from 
Silvaco  for  further  analysis. 

The  process  begins  with  the  DeckBuild  graphical  user 
interface  and  the  creation  of  an  input  file.  ATLAS  is  the 
module  used  in  this  research  and  the  module  is  adequate  and 
effective  for  all  of  the  work  that  was  accomplished .  The 
physical  parameters  and  structural  dimensions  were  set  for 
the  input  file;  to  be  used  globally  throughout  the  model, 
in  conjunction  with  a  doping  level  and  MODEL  command 
setting.  In  the  MODEL  command,  the  lattice  temperature  is 
set  for  the  ATLAS  run.  At  this  point,  the  program  is  ready 
to  enter  the  recording  stage  of  the  process. 

The  PROBE  function  was  used  extensively  to  gather  the 
data  required  for  proper  analysis.  This  function  was  used 
to  record  the  temperature,  the  hole  concentration,  and  the 
electron  concentration.  An  iterative  stepping  of  minimal 
bias  is  applied  to  allow  for  the  software  to  calculate  the 
properties  of  the  material.  The  STRUCTURE  (.STR)  file  is 
saved  and  the  probed  values  are  logged  into  a  LOG  (.log) 
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file.  At  this  point,  the  program  is  ready  to  quit  the 
ATLAS  run  at  the  given  temperature  and  begin  stepping  on  to 
the  next  temperature  to  begin  the  process  anew.  The  new 
temperature  is  set  in  the  MODEL  command.  The  LOG  file  is 
then  appended  to  record  the  new  data  from  the  updated 
temperature.  The  process  repeats  itself  until  the  QUIT 
command  is  called  to  end  the  process.  This  raw  data  can 
then  be  exported  into  Microsoft  Excel  or  MATLAB  to  be 
further  categorized  and  analyzed. 

B.  DECKBUILD  ABBREVIATED  SOURCE  CODE 

go  atlas 

set  dopant  =  boron 

#  Filneame 

set  datadump  =  ThesisSl 

set  e-alpha  =  0.00047 

set  e-beta  =  636 

set  e-band  =  1.42 

set  Eae  =  0.045 

set  Pae  =  0.02 

mesh  space. mult=l 

x.mesh  loc  =  0.00  spac  =  0.10 

x. mesh  loc  =  1  spac  =  0.10 

y. mesh  loc  =  0.00  spac  =  0.02 
y.mesh  loc  =  1  spac  =  0.02 

#  placing  anodes  and  cathodes 
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region  number  =  1  material  =  Si 

electrode  name  =  anode  top 

electrode  name  =  cathode  bottom 


material  region=l  egalpha=$ "e-alpha"  egbeta=$ "e-beta" 
EG300=$ "e-band"  EAB  =$"Eae"  EDB  =$"Pae" 


doping  reg 
model  temp 


1  uniform  $"dopant"  conc=lel7 
100  incomplete  conmob  fermi 


PROBE  NAME  = 
PROBE  NAME  = 
y  =  0.2 
PROBE  NAME  = 
0.05  y  =  0.2 


Temperature_ (k)  LAT . TEMP  x  =  0.05  y  =  0.2 
Hole_Concentrat ion_ ( 1 /cmA3 )  P.CONC  x  =  0.05 

Electron_Concentration_(l/cmA3)  N.CONC  x  = 


solve  init 

solve  vanode  =  0.1 

solve  vanode  =  0.2 

save  outf  =  $ "datadump" . str 

save  outf  =  $ "datadump" 60 . str 

log  outfile  =  $ "datadump" . log 

solve  vanode  =  1 


go  atlas 

mesh  inf ile=$ "datadump" . str 

material  region=l  egalpha=$ "e-alpha"  egbeta=$ "e-beta" 
EG300=$ "e-band"  EAB  =$"Eae"  EDB  =$"Pae" 

models  temp  =  110  incomplete  conmob  fermi 

PROBE  NAME=  Temperature. ( k )  LAT . TEMP  x=0.05  y=0 . 2 
PROBE  NAME=  Hole_Concent rat ion_ ( 1 / cmA 3 )  P.CONC  x=0.05 
y=0 . 2 

PROBE  NAME  =  Elect ron.Concent rat ion_ ( 1 / cmA 3 )  N.CONC  x  = 
0.05  y  =  0.2 

solve  init 

save  out f=$ "datadump" . str 

solve  vanode  =  0.1 

solve  vanode  =  0.2 

save  outf  =  $ "datadump" 1 1 0 . str 
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log  out f ile=$ "datadump" . log  append 
solve  vanode=l 
go  atlas 

mesh  inf ile=$ "datadump" . str 

material  region=l  egalpha=$ "e-alpha"  egbeta=$ "e-beta" 
EG300=$ "e-band"  EAB  =$"Eae"  EDB  =$"Pae" 

models  temp  =  120  incomplete  conmob  fermi 

PROBE  NAME=  Temperature. ( k )  LAT . TEMP  x=0.05  y=0 . 2 
PROBE  NAME=  Hole_Concent rat ion_ ( 1 / cmA 3 )  P . CONC  x=0.05 
y=0 . 2 

PROBE  NAME  =  Elect ron.Concent rat ion_ ( 1 / cmA 3 )  N.CONC  x 
0.05  y  =  0.2 

solve  init 

save  out f=$ "datadump" . str 

solve  vanode  =  0.1 

solve  vanode  =  0.2 

save  outf  =  $ "datadump" 12 0 . str 

log  out file=$ "datadump" . log  append 

solve  vanode=l 


go  atlas 

mesh  inf ile=$ "datadump" . str 

material  region=l  egalpha=$ "e-alpha"  egbeta=$ "e-beta" 
EG300=$"e-band"  EAB  =$"Eae"  EDB  =$"Pae" 


models  temp  =  1010  incomplete  conmob  fermi 

PROBE  NAME=  Temperature. ( k )  LAT . TEMP  x=0.05  y=0 . 2 
PROBE  NAME=  Hole_Concent rat ion_ ( 1 / cmA 3 )  P . CONC  x=0.05 
y=0 . 2 

PROBE  NAME  =  Elect ron.Concent rat ion_ ( 1 / cmA 3 )  N.CONC  x 
0.05  y  =  0.2 

solve  init 
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save  out f =$ "datadump" . str 

solve  vanode  =  0.1 

solve  vanode  =  0.2 

save  outf  =  $ "datadump" 1 0 1 0 . str 

log  out file=$ "datadump" . log  append 

solve  vanode=l 

go  atlas 

mesh  inf ile=$ "datadump" . str 

material  region=l  egalpha=$ "e-alpha"  egbeta=$ "e-beta" 
EG300=$"e-band"  EAB  =$"Eae"  EDB  =$"Pae" 


models  temp  =  1020  incomplete  conmob  fermi 

PROBE  NAME=  Temperature. ( k )  LAT . TEMP  x=0.05  y=0 . 2 
PROBE  NAME=  Hole_Concent rat ion_ ( 1 / cmA 3 )  P . CONC  x=0.05 
y=0 . 2 

PROBE  NAME  =  Electron_Concentration_(l/cmA3)  N.CONC  x 
0.05  y  =  0.2 

solve  init 

save  out f=$ "datadump" . str 

solve  vanode  =  0.1 

solve  vanode  =  0.2 

save  outf  =  $ "datadump" 1 02 0 . str 

log  out file=$ "datadump" . log  append 

solve  vanode=l 

go  atlas 

mesh  inf ile=$ "datadump" . str 

material  region=l  egalpha=$ "e-alpha"  egbeta=$ "e-beta" 
EG300=$"e-band"  EAB  =$"Eae"  EDB  =$"Pae" 


models  temp  =  1030  incomplete  conmob  fermi 

PROBE  NAME=  Temperature. ( k )  LAT . TEMP  x=0.05  y=0 . 2 
PROBE  NAME=  Hole_Concent rat ion_ ( 1 / cmA 3 )  P . CONC  x=0.05 
y=0 . 2 

PROBE  NAME  =  Elect ron.Concent rat ion_ ( 1 / cmA 3 )  N.CONC  x 
0.05  y  =  0.2 

solve  init 
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save  out f =$ "datadump" . str 

solve  vanode  =  0.1 

solve  vanode  =  0.2 

save  outf  =  $"datadump"1030 . str 

log  out file=$ "datadump" . log  append 

solve  vanode=l 

go  atlas 

mesh  inf ile=$ "datadump" . str 

material  region=l  egalpha=$ "e-alpha"  egbeta=$ "e-beta" 
EG300=$"e-band"  EAB  =$"Eae"  EDB  =$"Pae" 


models  temp  =  1040  incomplete  conmob  fermi 

PROBE  NAME=  Temperature. ( k )  LAT . TEMP  x=0.05  y=0 . 2 
PROBE  NAME=  Hole_Concent rat ion_ ( 1 / cmA 3 )  P . CONC  x=0.05 
y=0 . 2 

PROBE  NAME  =  Electron_Concentration_(l/cmA3)  N.CONC  x 
0.05  y  =  0.2 

solve  init 

save  out f=$ "datadump" . str 

solve  vanode  =  0.1 

solve  vanode  =  0.2 

save  outf  =  $"datadump"1040 . str 

log  out file=$ "datadump" . log  append 

solve  vanode=l 

go  atlas 

mesh  inf ile=$ "datadump" . str 

material  region=l  egalpha=$ "e-alpha"  egbeta=$ "e-beta" 
EG300=$"e-band"  EAB  =$"Eae"  EDB  =$"Pae" 


models  temp  =  1050  incomplete  conmob  fermi 

PROBE  NAME=  Temperature. ( k )  LAT . TEMP  x=0.05  y=0 . 2 
PROBE  NAME=  Hole_Concent rat ion_ ( 1 / cmA 3 )  P . CONC  x=0.05 
y=0 . 2 

PROBE  NAME  =  Elect ron.Concent rat ion_ ( 1 / cmA 3 )  N.CONC  x 
0.05  y  =  0.2 

solve  init 
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save  out f =$ "datadump" . str 

solve  vanode  =  0.1 

solve  vanode  =  0.2 

save  outf  =  $"datadump"1050 . str 

log  out file=$ "datadump" . log  append 

solve  vanode=l 

C.  PROCEDURAL  GUIDE  TO  CREATING  A  MOVIE 

The  ability  to  visualize  the  effects  of  temperature  on 
a  semiconductor  cell  is  a  very  impressive  method  for  illus¬ 
trating  the  concepts  involved.  For  this  purpose,  Silvaco 
has  the  capability  to  save  STRUCTURE  files  from  the 
DeckBuild  computations.  These  files  are  repositories  of 
information  for  the  semiconductor  cell  at  given  conditions. 
These  files  can  be  used  to  display  the  carrier  concentra¬ 
tion,  along  with  a  host  of  other  material  and  electrical 
properties.  To  access  these  files,  Silvaco  provides  a 
stand-alone  viewer  called  TONYPLOT . 

Multiple  STRUCTURE  files  can  be  loaded  into  TONYPLOT, 
and  the  display  drop  down  menu  item  can  be  used  to  display 
the  carrier  concentration  of  the  files.  All  of  these  files 
must  be  selected  because  the  programs  sets  the  displayed 
picture  relative  to  the  structure  windows  selected.  While 
there  is  a  effective  movie  animator  built  into  Silvaco,  the 
created  movie  is  not  exportable  into  other  media  format, 
thus  each  individual  file  must  be  screen  captured  into  a 
standard  graphical  file  type  such  as  (.jpg)  or  (.bmp),  for 
future  manipulation.  The  files  can  then  be  organized  in 
commercial  moviemaker  software  and  formatted  for  further 
animation.  (Microsoft  Office  PowerPoint  was  used  to  animate 
the  slides  into  a  movie.)  The  process  is  labor  intensive. 


53 


an 


illustrative  tutorial  of 


but  very  effective  in  creating 
the  desired  simulation  effect. 
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APPENDIX  C.  PAPER  PRESENTED  AND  ACCEPTED  AT  THE 
45™  IEEE  INTERNATIONAL  MIDWEST  SYMPOSIUM  ON 
CIRCUITS  AND  SYSTEMS  IN  TULSA,  OK  ON  AUG  4-7,  2002 


This  appendix  contains  the  "Comparison  of  Analytic  and 
Numerical  Models  with  Commercially  Available  Simulation 
Tools  for  the  Prediction  of  Semiconductor  Freeze-out 
Exhaustion"  paper  from  the  4-7  August  2002  45th  IEEE  Midwest 
Symposium  on  Circuits  and  System  at  Oklahoma  State 
University  in  Tulsa,  Oklahoma.  The  paper  was  a  condensed 
report  on  the  work  researched  in  this  thesis. 
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Comparison  of  Analytic  and  Numerical  Models  with 
Commercially  Available  Simulation  Tools  for  the  Prediction 
of  Semiconductor  Freeze-out  and  Exhaustion 

Ron  Pieper,  Sherif  Michael,  Derek  Reeves 
Naval  Postgraduate  School 
Department  of  Electrical  and  Computer  Engineering 

Monterey  CA  93943 
Email:  rjpieper@nps.navy.mil 


Abstract  -  Currently,  commercial 
software  packages,  such  as 
available  through  Silvaco 
International [1] ,  are  well 
designed  to  solve  the 
electron/hole  transport  problem. 
This  type  of  calculation  is 
usually  required  to  predict  the 
device  IV  characteristic. 
Surprisingly,  using  the  same 
package,  to  obtain  a  temperature 
dependent  plot  for  majority 
carrier  concentration  for  a 
uniform  semiconductor  requires  a 
somewhat  complicated  procedure 
[2]  .  Our  paper  will  present  an 
efficient  novel  way  of  obtaining 
this  curve  from  the  Silvaco 
International  software  and 
compare  the  results  with  a 
proposed  one  dimensional  single¬ 
equation  analytic  model  and  a 
numerical  model  that  predict  the 
temperature  dependence  for 
majority  concentration  in  all 
regimes  [3]  [4]  . 


Introduction 


Integrated  Circuits  are 
specified  to  operate  between 
designated  temperature  limits. 
The  circuit  designer  selects  the 
doping  level  or  levels  and 


typically  assumes  that  the 
dopants  are  approximately  100% 
ionized,  i.e.  exhaustion  and  the 
temperature  are  not  too  high. 
There  can  be  a  significant 
impact  on  the  value  of  a 
plethora  of  device  parameters 
such  as  depletions  widths  and/or 
FET  threshold  voltages  if  the 
assumption  is  violated.  In  the 
domain  that  the  temperature  is 
too  low  the  percentage 
ionization  of  dopant  or  dopants 
will  be  significantly  less  than 
100%.  The  value  for  the 
majority  carrier  concentration 
is  depressed  significantly  below 
the  design  value.  On  the  other 
hand  if  the  temperature  is  too 
high  the  thermal  generation 
effect  causes  the  majority 
carrier  concentration  to  become 
excessively  higher  than  the 
design  value  in  what  is  called 
the  intrinsic  temperature 
regime.  The  exhaustion  regime 
lies  between  this  two  extremes, 
intrinsic  and  freeze-out.  It  is 
well  known  that  for  a  multiple 
impurity  dopant  process  this 
"simple"  three-regime 
description  can  be  inadequate. 
What  is  important  to  the 
designer  is  the  plot  of  the 
majority  carrier  versus  the 
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versus 


temperature,  or  what  is  more 
commonly  done,  a  plot  of 
majority  carrier  concentration 
versus  reciprocal  of  the 
temperature . 

To-date  most  "analytic"  methods 
for  determining  dominant 

features  in  such  plots  make  use 
of  multiple  ad-hoc  arguments, 
which  taken  one  at  a  time, 
applies  in  only  two  of  the  three 
regimes  mentioned  above.  [5,  6, 

7]  Numerical  methods  based  on 
one-dimensional  analysis  are 
applicable  and  provide 

significant  flexibility  in  terms 
of  making  predictions  when  there 
are  multiple  dopants  and  when 
taking  into  account  2nd  order 
effects.  [4]  It  can  be  argued 
that  numerical  methods  are 
generally  going  to  be  applicable 
over  a  wider  range  of  problems 
than  what  can  be  solved  with 
exact  analysis.  Two  and  three 
dimensional  numerical  finite 
element  methods  provided  in  CAD 
based  commercial  packages  have 
been  widely  employed  by  the 
engineers  in  the  semiconductor 
processing  industry.  The  1-D 
results  both  numerical  and 
analytical  are  compared  with 
predictions  from  such  a  package 
of  tools  from  Silvaco  Int.  [1] 
See  Figure  (1)  .  For  details  on 
numerical  algorithm  please  see 
references  [3]  and  [4] . 

2 .  Basic  Analysis 

The  physical  basis  [1]  for 
analysis  assumes  that  there  is 
sufficient  spatial  uniformity  in 
the  doping  profiles  in  order  to 
apply  the  condition  of  local 
charge  neutrality. 


respectively,  while  Na  and  Nd 

are  the  ionized  acceptors  and 
donors  atom  concentrations 
respectively.  The  summations  run 
over  the  numbers  of  impurities. 
The  approach  presented  is 
facilitated  by  defining  the 
parameter  Z  where: 

Z  =  exp[(£R. -EF)/kT^  (2) 


where  k  is  Boltzman's  constant, 
Ef  is  the  Fermi  level  and  EFi  is 
the  intrinsic  Fermi  level.  In 
non-degenerate  cases,  i.e.  where 
the  Boltzman  approximation  can 
be  applied  to  simplify  the  exact 
expression  for  the  carrier 
concentrations,  the  Z  parameter 
is  directly  proportional  to  the 
equilibrium  hole  concentration. 
The  constant  of  proportionality 
is  the  temperature  dependent 
intrinsic  concentration.  It 

turns  out  that  Eq.  (1)  can  be 
expressed  as  P(Z) .  Numerical 
methods  based  on  successive 
substitution  in  solving  for  the 
zeros  in  expression  (1)  were 
found  to  be  unreliable  in 
producing  a  convergent  solution. 
However  the  method  proposed  here 
is  derived  from  the  numerical 
scheme  known  as  interval 
bisection  [8]  and  it  was  found 
work  well  with  variety  of 
combinations  of  profiles  and 
other  nonstandard  conditions 
The  effective  density  of  states, 
which  is  derived  in  most 
introductory  semiconductor 

texts,  is  given  by: 


N„ 


(2mnnkT) 
~h2 


3/2 


(3a) 


(  *  7 

3/2 

P(T)  =  Po~n0 

-I  A/ 

+2X 

=  0  (1) 

II 

* 

(3b) 

l  mnJ 

where  T 

is 

the 

Kelvin 

temperature. 

p0  and  nQ 

are  the 

where  m*  is 

the  effective 

mass 

equilibrium 

hole 

and  electron 

for  the  respective  bands 

[5]  . 

carrier 

concentrations 

In  order  to  facilitate  casting 
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Eq.  (1)  into  a  form  dependent  on 
the  Z-parameter  the  well-known 
relations  for  intrinsic 
concentration  and  Fermi  level 
are  used. 


n,  =  jKN, 


(4) 


Ec.  =  Ec  +  E"  +  kT  In 


JFi 


rN_ \  A 


(5) 


Silicon)  are  the  several 
occupation  degeneracy,  for 
acceptors  and  donors. 


respectively.  The  donor  and 
acceptor  concentrations  with  the 
over-bars  are  normalized  by  the 
intrinsic  concentration.  The 
energy  levels  for  the  impurities 
Ea,  Ed  are  defined  relative  to 
the  valence  band,  i.e., 

Ed=Ec~E'd  O) 


Following  the  standard  methods 
[5,  6,  7]  of  expressing  carrier 

concentrations  in  terms  of 
intrinsic  parameters  Eqs  (4-5), 
it  follows  that: 

n, 

no=—  (6a) 

0  Z 

Pa  =  niZ  (6b) 


where  E,  is  the  standard  cited 

value  for  donors  measured  with 
respect  to  the  conduction  band. 
This  analysis  will  also  be  the 
basis  for  an  analytic 

description  provided  in  the  next 
section . 

3.  Analytic  Solution, 
single  impurity 


Making  use  of  the  Fermi 
probability  distributions  [5]  to 
predict  ionization  levels  of  the 
donors  and  acceptors,  i.e., 

Na  and  N j  ,  and  substitution  of 

Eqs  (6)  and  Eq.  (1)  and  then 
dividing  by  the  intrinsic 
concentration  [3]  leads  to  a 
condition  on  P  (Z,  T)  : 


AT, 

P(Z)=0=Z+%  1 


+Kdr 


K  AT 

-Z'-T-Ea  ,7) 

1 7+k,z 


in  which  the  over¬ 
indicates  division 
the  temperature 

constants. 

-bar  notation 
by  na  ,  and 

dependent 

(EFI-Ed)/kT 

Kd,=gde 

(8a) 

al 

(8b) 

are  needed  to 

partial  ionization. 

characterize 
Also  Ma  and 

Md  are  the  numbers  of  impurity 
components,  and  ga  (4  typical  for 
Silicon)  and  gd  (2  typical  for 


Situation  is  described  by  the 
special  case  that  Ma=l  and  Md=0 
the  summation  subscript  in  Eq 
(7)  can  be  dropped  to  lighten 
the  notation.  Straightforward 
algebra  leads  to  a  cubic 
equation  in  Z. 


z3  +  —  z2  - 

K„ 


Na 


v1  +  ^y 


z- 


0  (10) 


The  corresponding  result  for  a 
single  impurity  N-type  is 
exactly  the  same  form,  producing 
a  cubic  in  Y=l/Z  with  revised 
coefficients  obtained  by  letting 
the  'a'  subscript  (for  acceptor) 
be  replaced  with  the  'd' 
subscript  (for  donor)  .  The 

solution  can  be  defined  in  terms 
of  coefficients  for  the 
reference  cubic  equation: 

x3  +  axx2  +  a2x  +  a3  =  0  (li) 


To  facilitate  representation  of 
the  solution  the  following 
intermediate  parameterization  of 
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the  problem  is  commonly  taken 
[9]  . 

Q  =  (3a2-af)/9  (12a) 

R  =  {9axa2  -  21  a2  -  2a13)/54  (12b) 

D  =  QijrR2  (12c) 


where,  D,  which  is  referred  to 
as  the  "discriminant"  for  the 
cubic  problem  will  dictate  the 
type  of  solutions  possible.  As 
per  the  fundamental  theorem  of 
algebra  there  will  be  three 
roots.  For  D  negative  all  roots 
are  real  while  for  D  positive 
only  one  root  is  real  while  the 
other  two  are  complex  conjugate. 
If  D  is  zero  there  will  be 
repeated  root.  For  the  problem 
being  evaluated  it  can  be  shown 
using  symbolic  mathematical 
methods  that  Q  and  D  will  always 
be  negative.  Furthermore  it  can 
then  be  shown  that  of  the  three 
roots  only  one  is  positive  and 
therefore  physically  acceptable. 
That  root  will  be  predicted  by 
the  following  recipe  [9] 


Z  =  x  = 


2^PQcos(OI3) 


(13) 


where  6  =  COS  1 1 R  /  yJ~Q3  j  •  Once 

Z  is  determined  the  majority 
carrier  concentration,  p0,  can  be 
predicted  from  Eq  (6b) . 


purposes  of  estimating 

incomplete  ionization. 

Due  to  lack  of  space  it  is  not 
possible  to  cover  all  the 
details  in  the  algorithm 
involving  Silvaco  Tools.  Here 
are  some  critical  points. 
First,  the  Silvaco  predictions 
are  obtained  from  running  the 
software  over  2  dimensional  grid 
covering  a  2 
uniform  and 

semiconductor.  A 
is  chosen  in  the 
grid  and  the  majority  carrier 
concentration  is  recorded.  This 
is  done  in  an  automated  fashion 
for  each  temperature  in  a  long 
list  of  temperatures  covering 
the  range  of  interest.  The 

Silvaco  based  recordings 

described  in  the  next  section 
are  obtained  via  this  procedure. 


of  points 
dimensional, 
equilibrium 
center  point 


5 .  Examples 

Figure  2,  which  has  on  it  6 
plots  of  loglO  (p0)  versus 
1000x(l/T)  for  various  conditions 
(with  and  without  band-gap 
narrowing  (BGN)  and  from  various 
sources  (numerical,  analytic, 
and  Silvaco) .  The  model  for 
band-gap  narrowing  [6]  used  here 
is  . 


(T)  =  ESO 


aT 2 
(T  +  J3) 


(14) 


4 .  Silvaco  Int .  based 
Simulation 

Typically,  semiconductor  device 
problems  of  interest  involve 
conditions  such  as  external 
voltages  and/or  spatial  non¬ 
uniformity  in  the  semiconductor. 
Under  these  more  general 
conditions,  a  solution  can  be 
obtained  by  solving  Poisson' s 
equation.  This  is  done  with  the 
Silvaco  tool  known  as  ATLAS  [2] . 
It  is  possible  to  incorporate 
Fermi-Dirac  statistics  for 


The  various  curves  are  coded  by 
number  with  the  description  in 
the  legend.  In  Table  1 
recommended  values  for  Silicon 
including  BGN  parameters  a  and  (3 
are  provided.  Note,  as 
expected,  curves  from  the 
numerical  and  analytic  methods 
are  indistinguishable. 
Furthermore,  as  expected  the 
impact  of  BGN  is  apparent  in  the 
curves  at  the  higher 
temperatures.  The  difference 
between  Silvaco  predictions  and 
the  1-dimensional  modeling 
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efforts  are  slight  but  noted  for 
future  investigation. 


Parameter 

Value 

Ea 

0.045  ev 

Na 

1017  #/cm3 

Ego 

1 . 1 lev 
(Silicon) 

ga 

4 

BGN  a  [6] 

0 . 00047 

BGN  P  [6] 

636.0 

Table  1 :  Parameters  used  in 

Simulations 

6.  Conclusions 


Both  analytic  one-dimensional 


models  and  numerical  one¬ 
dimensional  algorithms  are 
essentially  indistinguishable  in 
terms  of  predictions.  As 
expected,  when  BGN  effects  are 
included  in  the  models  the 
impact  was  observed  at  the 
higher  temperatures.  The 
difference  between  Silvaco 
predictions  based  two- 
dimensional  modeling  and  the 
one-dimensional  solutions  were 


observed  to  be  slight  but  worthy 
of  future  investigation. 


FreffiBoi  exhaustion  regrrespreddions 


Figure  1:  Testing  Flowchart 


0.500  2.500  4.500  6.500  8.500  10.500 

1000  X  1/T 


Figure  2  :  Composite  Plot  of 

Results 
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